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We demonstrate the accuracy of the hypernetted chain closure and of the mean-field approxi- 
mation for the calculation of the fluid-state properties of systems interacting by means of bounded 
and positive-definite pair potentials with oscillating Fourier transforms. Subsequently, we prove the 
validity of a bilinear, random-phase density functional for arbitrary inhomogeneous phases of the 
same systems. On the basis of this functional, we calculate analytically the freezing parameters 
of the latter. We demonstrate explicitly that the stable crystals feature a lattice constant that 
is independent of density and whose value is dictated by the position of the negative minimum 
of the Fourier transform of the pair potential. This property is equivalent with the existence of 
clusters, whose population scales proportionally to the density. We establish that regardless of the 
form of the interaction potential and of the location on the freezing line, all cluster crystals have 
a universal Lindemann ratio Lf = 0.189 at freezing. We further make an explicit link between the 
aforementioned density functional and the harmonic theory of crystals. This allows us to establish 
an equivalence between the emergence of clusters and the existence of negative Fourier components 
of the interaction potential. Finally, we make a connection between the class of models at hand and 
the system of infinite-dimensional hard spheres, when the limits of interaction steepness and space 
dimension are both taken to infinity in a particularly described fashion. 
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I. INTRODUCTION 

Cluster formation in complex fluids is a topic that has 
attracted considerable attention recentlyjii^i^i^i^i^ii^'^ 
The general belief is that a short-range attraction in the 
pair interaction potential is necessary to initiate aggre- 
gation and a long-range repulsive tail in order to limit 
cluster growth and prevent phase separation. An al- 
ternative scenario for cluster formation pertains to sys- 
tems whose constituent particles interact by means of 
purely repulsive potentials. Cluster formation in this 
case is counterintuitive at first sight: why should par- 
ticles form clusters if there is no attraction acting be- 
tween them? The answer lies in an additional property 
of the effective repulsion, namely that of being bounded, 
thus allowing full particle overlaps. Though surprising 
and seemingly unphysical at first, bounded interactions 
are fully legitimate and natural as effective potentials'^ 
between polymeric macromolecular aggregates of low in- 
ternal monomer concentration, such as polymerS j^^'^^i^^ 
dendrimerSfi^ii^ microgelsr^iiiii^ or coarse-grained block 
copolymersii2i2£ The growing interest in this type of 
effective interactions is also underlined by the recent 
mathematical proof of the existence of crystalline ground 
states for such systemsi^i^^ 

Cluster formation in the fluid and in the crystal phases 
was explicitly seen in the system of penetrable spheres, 
following early simulation results^ and subsequent cell- 



model calculations.^'' Cluster formation was attributed 
there to the tendency of particles to create free space by 
forming full overlaps. The conditions under which ultra- 
soft and purely repulsive particles form clusters have been 
conjectured a few years ago^ and explicitly confirmed by 
computer simulation very recentlyi^ The key lies in the 
behavior of the Fourier transform of the effective interac- 
tion potential: for clusters to form, it must contain nega- 
tive parts, forming thus the class of Q^-interactions. The 
complementary class of potentials with purely nonnega- 
tive Fourier transforms, Q^, does not lead to clustering 
but to remelting at high densities i^^'^^i^^'^°i'^^'^^ An in- 
triguing feature of the crystals formed by Q^-systems is 
the independence of the lattice constant on density i^^i^^ 
a feature that reflects the flexibility of soft matter sys- 
tems in achieving forms of self-organization unknown to 
atomic onesi^^iM The same characteristic has recently 
been seen also in slightly modified models that contain 
a short-range hard corei^ In this work, we provide an 
analytical solution of the crystallization problem and of 
the properties of the ensuing solids within the frame- 
work of an accurate density functional approach. We 
explicitly demonstrate the persistence of a single length 
scale at all densities and for all members of the Q*-class 
and offer thus broad physical insight into the mecha- 
nisms driving the stability of the clustered crystals. We 
further establish some universal structural properties of 
all Q^-systems both in the fluid and in the solid state, 
justifying the use of the mean-field density functional 
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on which this work rests. We make a connection be- 
tween our results and the harmonic theory of solids in 
the Einstein-approximation. Finally, we establish a con- 
nection with suitably-defined infinite dimensional models 
of hard spheres. 

The rest of this paper is organized as follows: in Sec. 
ITTl we derive an accurate density functional by starting 
with the uniform phase and establishing the behavior of 
the direct correlation functions of the fluid with density 
and temperature. Based on this density functional, we 
perform an analytical calculation of the freezing charac- 
teristics of the Q^-systems by employing a weak approx- 
imation in Sec. IIIII The accuracy of this approximation 
is successfully tested against full numerical minimization 
of the functional in Sec. II VI In Sec.|V]the equivalence be- 
tween the density functional and the theory of harmonic 
crystals is demonstrated, whereas in Sec. I VII a connec- 
tion is made with inverse-power potentials. Finally, in 
Sec. IVIl] we summarize and draw our conclusions. Some 
intermediate, technical results that would interrupt the 
flow of the text are relegated in the Appendix. 



II. DENSITY FUNCTIONAL THEORY 

A. Definition of the model 

In this work, we focus our interest on systems 
of spherosymmetric particles interacting by means of 
bounded pair interactions v{r) of the form: 

v{r) = e(l){r/a), (1) 

with an energy scale e and a length scale a, and which 
fulfill the Ruelle conditions for stability!^ In Eq. (P) 
above, (/^{x) is some dimensionless function of a dimen- 
sionless variable and r denotes the distance between the 
spherosymmetric particles. In the context of soft matter 
physics, v{r) is an effective potential between, e.g., the 
centers of mass of macromolecular entities, such as poly- 
mer chains or dendrimers.^^ As the centers of mass of the 
aggregates can fully overlap without this incurring an in- 
finitely prohibitive cost in (free) energy, the condition of 
boundedness is fulfilled: 

< v{r) < K, (2) 

with some constant K. 

The interaction range is set by ct, typically the physical 
size (e.g., the gyration radius) of the macromolecular ag- 
gregates that feature v{r) as their effective interaction. In 
addition to being bounded, the second requirement to be 
fulfilled by the function (f){x) is that it decay sufficiently 
fast to zero as a; — > oo, so that its Fourier transform (f>{y) 
exists. In three spatial dimensions, we have 

4>{y) = An r dx ^^^^x^cl>{x) (3) 
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and, correspondingly, 

v{k) = ea^4>{ka) (4) 

for the Fourier transform v{k) of the potential, evaluated 
at wavenumber k. Our focus in this work is on systems 
for which v(k) is oscillatory, i.e., v{r) features positive 
and negative Fourier components, classifying it as a Q^- 
potentiaL^-£ Though this work is general, for the purposes 
of demonstration of our results, we consider a particular 
realization of Q*-potentials, namely the generalized ex- 
ponential model of exponent m, (GEM-m): 

v{r)^eexp[-{r/ar], (5) 

with TO = 4. It can be shown that all members of the 
GEM-TO family with to > 2 belong to the Q*-class, see 
Appendix A. 

A short account of the freezing and clustering behav- 
ior of the GEM-4 model has been recently publishedi^ 
Another prominent member of the family is the m = oo 
model, which corresponds to penetrable spheres^^i^ii^i 
with a finite overlap energy penalty e. Indeed, the explic- 
itly calculated phase behavior of these two show strong 
resemblances, with the phase diagram of both being dom- 
inated by the phenomenon of formation of clusters of 
overlapping particles and the subsequent ordering of the 
same in periodic crystalline arrangementsj^ii?''' In this 
work, we provide a generic, accurate, and analytically 
tractable theory of inhomogeneous phases of Q^-systems. 

B. The uniform fluid 

Let us start from the simpler system of a homogeneous 
fluid, consisting of N spherosymmetric particles enclosed 
in a macroscopic volume V . The structure and thermo- 
dynamics of the system are determined by the density 
p = N/V and the absolute temperature T or, better, 
their dimensionless counterparts: 

p* ^ pa^ 

and 

5 

e 

with Boltzmann's constant fee. As usual, we also in- 
troduce for future convenience the inverse temperature 
(3 = (/cbT)^^. We seek for appropriate and accurate clo- 
sures to the Ornstein-Zernike relation^ 

h{r) = c(r) +p j d-V'c(|r - r'\)h{r'), (6) 

connecting the total correlation function h(r) to the di- 
rect correlation function c(r) of the uniform fluid. One 
possibility is offered by the hypernetted chain closure 
(HNC) that reads as^ 

h{r) = exp [-I3v{r) + h{r) - c(r)] - 1. (7) 
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An additional closure, the mean-field approximation 
(MFA), was also employed and will be discussed later. 

Our solution by means of approximate closures was 
accompanied by extensive NVT-Monte Carlo (MC) sim- 
ulations. We measured the radial distribution func- 
tion g{r) = h{r) + 1 as well as the structure factor 
S{k) — 1 -f ph{k), where h{k) is the Fourier transform 
of h{r), to provide an assessment of the accuracy of the 
approximate theories. We typically simulated ensembles 
of up to 3 000 particles for a total of 150 000 Monte Carlo 
steps. Measurements were taken in every tenth step after 
equilibration. 

In Fig. [1] we show comparisons for the function g{r) 
as obtained from the MC simulations and from the HNC 
closure for a variety of temperatures and densities. It 
can be seen that agreement between the two is obtained, 
to a degree of quality that is excellent. Tiny deviations 
between the HNC and MC results appear only at the 
highest density and for a small region around r = for 
low temperatures, T* — 0.5. Otherwise, the system at 
hand is described by the HNC with an extremely high 
accuracy and for a very broad range of temperatures and 
densities. We note that, although in Fig. [1] we restrict 
ourselves to temperatures T* < 2.0, the quality of the 
HNC remains unaffected also at higher temperaturesi^ 

In attempting to gain some insight into the remark- 
able ability of the HNC to describe the fluid structure 
at such a high level of accuracy, it is useful to recast 
this closure in density-functional language. Following 
the famous Percus idea,^^'^^-'^'^''^^''^'^ the quantity pg{r) 
can be identified with the nonuniform density p{r) of an 
inhomogeneous fluid that results when a single particle 
is kept fixed at the origin, exerting an 'external' poten- 
tial Vcxtir) — v{r) on the rest of the system. Following 
standard procedures from density functional theory, we 
find that the sought-for density profile p{r) is given by 

p{r) = A-3 exp |/3m - (3v{r) - ^^^j ^ («) 

where A is the thermal de Broglie wavelength and p, the 
chemical potential associated with average density p and 
temperature T. Moreover, -Fox[p] is the intrinsic excess 
free energy, a unique functional of the density p{r). As 
such, -Fexl/o] can be expanded in a functional Taylor series 
around its value for a uniform liquid of some (arbitrary) 
reference density po- For the problem at hand, po = p is 
a natural choice and we obtain'*'^ 

/3i^ex[p] ^PFM - Y.-\ I I I d'^idV2 . . . d3r„ 
„=i ^- J J J 

X c^"^(ri,r2, . . . ,r„;p) 

X Ap(ri)Ap(r2)...Ap(r„), (9) 

where F^xip) denotes the excess free energy of the homo- 
geneous fluid, as opposed to that of the inhomogeneous 
fluid, Fcx[p], and 
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FIG. 1: Radial distribution functions g{r) of the GEM-4 
model as obtained by Monte Carlo simulation (points), the 
HNC-closure (solid lines) and the MFA (dashed lines), at var- 
ious temperatures and densities. For clarity, the curves on 
every panel have been shifted upwards by certain amounts, 
which are given below in square brackets, following the value 
indicating the density p*. (a) T* — 0.5 and densities, from 
bottom to top: p* = 0.5 [0]; p* = 1.0 [0.2]; p* = 1.5 [0.4]; 
p* = 2.0 [0.6]; p* = 2.5 [0.8]. (b) T* = 1.0 and densities, from 
bottom to top: p* = 1.0 [0]; p* = 2.0 [0.2]; p* = 3.0 [0.4]; 
p* = 4.0 [0.6]; p* = 5.0 [0.8]. (c) T* = 2.0 and densities, from 
bottom to top: p* = 2.0 [0]; p* = 4.0 [0.2]; p* = 6.0 [0.4]; 
p* = 8.0 [0.6]; p* = 10.0 [0.8]. 



Ap(x) = p(x) - p. 



(10) 
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The basis of the expansion given by Eq. ([9]) above is 
the fact that i^ex[p] is the generating functional for the 



hierarchy of the direct correlation functions (dcf 's) 



J") 



In particular, — Cq is the n-th functional derivative 
of the excess free energy with respect to the density field, 
evaluated at the uniform density p-M>^ 



(ri,r2. 



p(-}=p 
(11) 

As the functional derivatives are evaluated for a uni- 
form system, translational and rotational invariance re- 
duce the number of variables on which the n-th or- 



der dcf's Cq"'' depend; in fact, CQ^^(ri;p) is a position- 
independent constant and equals — /J/iox, where /icx is 
the excess chemical potential4i Similarly, Cp^^ (ri, r2; p) 
is simply the Ornstein-Zernike direct correlation func- 
tion c(|ri — r2|) entering in Eqs. ^ and ([7]) above. The 
HNC closure is equivalent to jointly solving Eqs. © 
and ([8]) by employing an approximate excess free energy 
functional i^exlp], arising by a truncation of the expan- 
sion of Eq. ([5]) at n = 2, i.e., discarding all terms with 
n > 3; this is the famous Ramakrishnan-Yussouff second- 
order approximation4Si^ Indeed, the so-called bridge 
function b(r) can be written as an expansion over in- 
tegrals involving as kernels all the Cg"-* with n > 3 and 
the HNC amounts to setting the bridge function equal to 
zero^^iil 

Whereas in many cases, such as the one-component 
plasma j^^i**^ and other systems with long-range interac- 
tions, the HNC is simply an adequate or, at best, a very 
good approximation, for the case at hand the degree of 
agreement between the HNC and simulation is indeed ex- 
tremely high. What is particularly important is that the 
accuracy of the HNC persists for a very wide range of 
densities, at all temperatures considered. This fact has 
far-reaching consequences, because it means that the cor- 
responding profiles /o(x) that enter the multiple integrals 
in Eq. ^ vary enormously depending on the uniform 
density considered. Thus, it is tempting to conjecture 
that for the systems under consideration (soft, penetra- 
ble particles at T* > 1), not simply the integrals with 
n > 3 vanish but rather the kernels themselves. In other 
words. 



(1) 



(ri,r2, . . . ,r„;p) 0, 



{n > 3). 



(12) 



The behavior of the higher-order dcf's is related to 
the density-derivative of lower-order ones through certain 
sum rules, to be discussed below. Hence, it is pertinent 
to examine the density dependence of the dcf c(r) of the 
HNC. In Fig. [5] we show the difference between the dcf 
chnc(^) and the mean-field approximation (MFA) to the 
same quantity: 



CMFA(f) = ~l3v{r). 



(13) 



Eq. (|13p above is meaningless if the pair potential v{r) 
diverges as r ^ 0, because c(r) has to remain finite at 
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FIG. 2: The main plots show the difference between the di- 
rect correlation function c(r) calculated in the HNC and its 
MFA- approximation, c(r) — —f3v{r), for a GEM-4 model, at 
various densities indicated in the legends. The insets show the 
MFA approximation for the same quantity, which is density- 
independent. Each panel corresponds to a different temper- 
ature: (a) T* = 0.5; (b) T* = 1.0; (c) T* = 2.0. These 
are exactly the same parameter combinations as the ones for 
which g{r) is shown in Fig. [T] 
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all r, as follows from exact diagrammatic expansions of 
the samcj^ In fact, the form c(r) ~ —/3v{r) denotes the 
large-r asymptotic behavior of c(r). In our case, how- 
ever, where v{r) lacks a hard core, the MFA-form for 
c(r) cannot be a priori rejected on fundamental grounds; 
the quantity —(3v(r) remains bounded as r — ^ 0. In fact, 
as can be seen in Fig. [21 the deviations between the MFA 
and the HNC-closure are very small for T* > 1. In addi- 
tion, the differences between chnc('") and Cmfa('') drop, 
both in absolute and in relative terms, as temperature 
grows, see the trend in Figs. [2{a)-{c). At fixed tempera- 
ture, the evolution of the difference with density is non- 
monotonic: it first drops as density grows and then it 
starts growing again at the highest densities shown in 
the three panels of Fig. [21 

Motivated by these findings, we employ now a second 
closure, namely the above-mentioned MFA, Eq. p^ . In- 
troducing the latter into the Ornstcin-Zcrnike relation, 
Eq. ([5]), we obtain the MFA-results for the radial distri- 
bution function g{r) that are also shown in Fig. [1] with 
dashed lines. Before proceeding to a critical comparison 
between the g{r) obtained from the two closures, it is 
useful to make a clear connection between the MFA and 
the HNC. 

As mentioned above, the members of the sequence 
of the n-th order direct correlation functions are not 
independent from one another; rather, they are con- 
strained to satisfy a corresponding hierarchy of sum rules, 
namely^ii^^i^ 



j3 ("+!)/ \ 

d'^rfcc;, '(ri, . . . ,rfc_i,rfe,rfc+i, . . . ,r„+i;p) = 

9c[)"^(ri, . . . ,rfc-i,rfc+i, . . . ,r„+i;p) 
dp 

In particular, for n = 2, we have 

J dV4^)(r,rMr-r'|;p) = 



(14) 



dV4')(/,r,|r'-r|;p) 

dc{r; p) 
dp ' 



(15) 



where we have used the translational and rotational in- 
variance of the fluid phase to reduce the number of argu- 
ments of the dcf's and we show explicitly the generic 
dependence of c(r) on p. In the MFA, one assumes 
c(r) — —f3v{r), with the immediate consequence 



dc{r) 



0. 



(16) 



Eqs. and imply that the integral of c"q^ with 
respect to any of its arguments must vanish for arbitrary 
density. As Cg has a complex dependence on its ar- 
guments, this is a strong indication that c"^ itself van- 
ishes. In fact, both for the Barrat-Hansen-Pastore factor- 
ization approximation for this quantit y^*^'^^ and for the 



alternative, Denton- Ashcroft fc-space factorization of the 
same,^^ the vanishing of the right-hand side of Eq. p5|) 



implies that Cp — 0. Now, if c"q' vanishes, so does also 
its density derivative and use of sum rule (|14p for n = 3 

implies c"q^ = 0. Successive use of the same for higher 
n-values leads then to the conclusion that in the MFA: 



.(3) 



„(") 



(ri,r2, . 



;p)-0, (n>3). 



(17) 



We can now see that the accuracy of the HNC stems 
from the fact that for these systems we can write 



c(r) = — /3w(r) -f £(r; p). 



(18) 



where e{r\ p) is a small function at all densities p, offering 
concomitantly a very small, and the only, contribution 
to the quantity dc(r) / dp. This implies that c"q^ itself is 
negligible by means of Eq. (fTSj) . Repeated use of Eq. (fH|) 
leads then to Eqs. (fT2l) and shows that the contributions 
from the n > 2 terms, that are ignored in the HNC, 
are indeed negligible. The HNC is, thus, very accurate, 
due to the strong mean-field character of the fiuids at 
hand." The deviations between the HNC and the MFA 
come through the function £(r; p) above. 

Let us now return to the discussion of the results for 
g{r) and c(r) and the relative quality of the two closures 
at various thermodynamic points. Referring first to Fig. 
mc), we see that at T* = 2.0 both the HNC and the 
MFA perform equally well. The agreement between the 
two (and between the MFA and MC) worsens somewhat 
at T* = 1.0 and even more at T* ^ 0.5. The MFA is, 
thus, an approximation valid for T* > 1, in agreement 
with previous resultsi^^ The reason lies in the accuracy 
of the low-density limit of the MFA. In general, as p ^ 0, 
c(r) tends to the Mayer function /(r) = exp[— /3u(r)] — 1. 
If 7"* ^ li one may expand the exponential to linear 
order and obtain /(r) = —(3v{r), so that the MFA can 
be fulfilled. 

In the HNC, it is implicitly assumed that all dcf's with 
n > 3 vanish. In the MFA, this is also the case. The 
two closures differ in one important point, though: in 
the HNC, the second-order direct correlation function is 
not prescribed but rather determined, so that both the 
'test-particle equation', Eq. ([S]), and the Ornstein-Zernike 
relation, Eq. ([6]), are fulfilled. In the MFA, it is a pri- 
ori assumed that c(r) — —(3v{r), which is introduced 
into the Ornstein-Zernike relation and thus h{r) is deter- 
mined. This is one particular way of obtaining g{r) in the 
MFA, called the Ornstein-Zernike route. Alternatively, 
one could follow the test-particle route in solving Eq. ^ 
in conjunction with Eq. ^ and the MFA-approximation, 
Eq. (fT5)) . In this case, the resulting expression for the to- 
tal correlation hmction h{r) in the MFA reads as 



h{r) — exp[—f3v{r) — f3p{h * v){r)] — 1, 



(19) 



with * denoting the convolution. Previous studies with 
ultrasoft systems have shown that the test-particle h{r) 
from the MFA is closer to the HNC-result than the MFA 
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result obtained from the Ornstein-Zernike router^i The 
discrepancy between the two is a measure of the approx- 
imate character of the MFA; were the theory to be ex- 
act, aU routes would give the same result. As a way to 
quantify the approximations involved in the MFA, let us 
attempt to impose consistency between the test-particle 
and Ornstein-Zernike routes. Since c(r) = —(3v(r) in 
this closure, the exponent in Eq. above is just the 
right-hand side of the Ornstein-Zernike relation, Eq. ([6]). 
Thus, if we insist that the latter is fulfilled, we obtain the 
constraint 



h{r) ~ exp[h{r)] — 1, 



(20) 



which is strictly satisfied only for h{r) = 0. However, 
as long as \h{r)\ < 1, one can linearize the exponential 
and an identity follows; the internal inconsistency of the 
MFA is of quadratic order in h{r) and it follows that 
the MFA provides an accurate closure for the systems at 
hand, as long as \h{r)\ remains small. This explains the 
deviations between MFA and MC seen at small r for the 
highest density at T* — 2.0, Fig. [T{c). The same eSect 
can also be seen in Fig. ^c) as a growth of the discrep- 
ancy between chnc(?') and cmfa(?') at small r- values for 
the highest density shown. In absolute terms, however, 
this discrepancy remains very small. Note also that dis- 
crepancies at small r-values become strongly suppressed 
upon taking a Fourier transform, due to the additional 
geometrical r^-factor involved in the three-dimensional 
integration. 

It can therefore be seen that the MFA and the HNC 
are closely related to one another: the HNC is so suc- 
cessful due to the strong mean-field character of the sys- 
tems under consideration. This fact has also been estab- 
lished and extensively discussed for the case of the Gaus- 
sian modelp^i^ i.e., the m = 2 member of the GEM-to 
class. Once more, the HNC and the MFA there are very 
accurate for high densities and/or temperatures, where 
h{r) = and the system's behavior develops similari- 
ties with an 'incompressible ideal gas',~ in full agree- 
ment with the remarks presented above. Subsequently, 
the MFA and HNC closures have been also successfully 
applied to the study of structure and thermodynamics of 
binary soft mixtures i^'^'^'S^i^'^ 

A crucial difference between the Gaussian model, 
which belongs to the (5~'"-class, and members of the Q^- 
class, which are the subject of the present work, lies in the 
consequences of the MFA-closure on the structure factor 
S{k) of the system. Since c(r) = —(3v{r), the Ornstein- 
Zernike relation leads to the expression 



Sik) = 



1 



1 + f3pv{k) ' 



(21) 



Whereas for Q+-potentials S{k) is devoid of pronounced 
peaks that exceed the asymptotic value S{k —f oo) = 1, 
for Q^-systems a local maximum of S{k) appears at the 
value fc, for which v{k) attains its negative minimum. 
In Fig. [3] we show representative results for the system 
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FIG. 3: The structure factor S{k) of the GEM-4 model at 
temperature T* = 1.0 and various densities, indicated be- 
low. For clarity, the curves have been shifted vertically by 
amounts shown in the square brackets, following the numbers 
that indicate the values of the density p*. From bottom to 
top: p* = 1.0 [0]; p* = 2.0 [0.5]; p* = 3.0 [1.0]; p* = 4.0 [1.5]; 
p* = 5.0 [2.0]. The points are results from Monte Carlo sim- 
ulations and the dashed lines from the HNC. As the HNC- 
and MFA-curves run very close to each other, we show the 
MFA-result by the solid curve only for the highest density, 
p* = 5.0. 



at hand, where it can also be seen that the HNC and 
the MFA yield practically indistinguishable results. In 
full agreement with the MC simulations, the location of 
the main peak of S{k) is density-independent, a feature 
unknown for usual fluids, having its origin in the fact 
that c(r) itself is density independent. 

Associated with this is the development of a A- 
line^^i^i^'^ also known as Kirkwood instability,^^ on 



which the denominator in Eq. (j2ip vanishes at /c* and 
thus S'(/c*) oo. The locus of points {px,T\) on 
the density-temperature plane for which this divergence 
takes place is, evidently, given by 



-4i{kf,cr). 



(22) 



In the region p> p\ (equivalently: T < T\) on the {p, T)- 
plane, the MFA predicts that the fluid is absolutely un- 
stable, since the structure factor there has multiple diver- 
gences and also develops negative parts. This holds only 
for Q^-systems; for Q^-ones the very same line of argu- 
mentation leads to the opposite conclusion, namely that 
the fluid is the phase of stability at high densities and/or 
temperatures. The latter conclusion has already been 
reached by Stillinger and coworker a^"^'^^i^^'^^i^^ in their 
pioneering work of the Gaussian model in the mid-1970's, 
and explicitly confirmed by extensive theory and com- 
puter simulations many years later i^^i^^i^^i^'^'^-'^i'^^ How- 
ever, Stillinger's original argument was based on duality 
relations that are strictly fulfilled only for the Gaussian 
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model, whereas the MFA-arguments are quite general. 



C. Nonuniform fluids 

Having established the validity of the MFA for vast 
domains in the phase diagram of the systems under con- 
sideration as far as the uniform fluids are concerned, we 
now turn our attention to nonuniform ones. Apart from 
an obvious general interest in the properties of nonuni- 
form fluids, the necessity to consider deviations from ho- 
mogeneity in the density for Q^-models is dictated by 
the A-instability mentioned above: the theory of the uni- 
form fluid contains its own breakdown, thus the system 
has to undergo a phase transformation to a phase with a 
spontaneously broken translational symmetry. Whether 
this transformation takes place exactly on the instabil- 
ity line or already at densities p < p\ (or temperatures 
T > T\) and which is the stable phase are some of the 
questions that have to be addressed. Density-functional 
theory of inhomogeneous systems is the appropriate the- 
oretical tool in this direction. 

Let us consider a path Px{^) in the space of density 
functions, which is characterized by a single parameter 
X; this path starts at some reference density Pr(r) and 
terminates at another density p(r). The uniqueness of 
the excess free energy functional and its dependence on 
the inhomogeneous density field p{y) allow us to inte- 
grate dFc-x[p\/ dx along this path, obtaining Fcx[/o], pro- 
vided that -Fox[Pr] is known. A convenient paramctriza- 
tion reads as Pxi'^) — Pr(r) + x[p(r) — Pr(r)], with x = 
corresponding to Pr(r) and x = 1 to p(r). The excess 
free energy of the final state can be expressed as^ 



l3F,^[p] = l3F,^[p,] 



dx / dVc(i)(r;[p^])Ap(r), (23) 



where c^^^(r; [p^]) denotes the first functional derivative 
of the quantity — /3i^cx[p] evaluated at the inhomoge- 
neous density Px(r), and Ap(r) — p{r) — Pr{r). Since 
c^^^{r; [px]) is in its own turn a unique functional of the 
density profile, repeated use of the same argument leads 
to a functional Taylor expansion of the excess free energy 
around that of a reference system, an expansion that ex- 
tends to infinite order. For the particular choice of a 
uniform reference system, /9r(r) = p, we obtain then the 
Taylor series of Eq. In general, however, the reference 
system does not have to have the same average density 
as the final one, hence the uniform density p there must 
be replaced by a more general quantity pQ. 

The usefulness of Eq. ^ in calculating the free en- 
ergies of extremely nonuniform phases, such as crystals, 
is limited both on principal and on practical grounds. 
Fundamentally, there is no small parameter guiding 
such an expansion, since the differences between the 
nonuniform density of a crystal and that of a fiuid are 
enormous; the former has extreme variations between 



lattice- and interstitial regions. Hence, the very con- 
vergence of the series is in doubt. In practice, the di- 
rect correlation functions for n — 3 are very cumber- 
some to calculate- - and those for n > 4 are practically 
unknown."*^ The solution is either to arbitrarily terminate 
the series at second order— or to seek for nonperturbative 
functionalS ]^^'''"°i'^"'^''^^i'^^''^"^i'^^ In our case, however, things 
are different because, for the systems we consider, we 
have given evidence that the dcf's of order n > 2 are 
extremely small and we take them at this point as van- 
ishing. Then, the functional Taylor expansion of the free 
energy F^x [p] terminates (to the extent that the approxi- 
mation holds) at second order. The Taylor series becomes 
a finite sum and convergence is not an issue any more. 

Let us, accordingly, expand -Fcx[p] around an arbi- 
trary, homogeneous reference fluid of density po — Nq/V, 
taking into account that the volume V is flxed but the 
system with density p(r) contains N particles, whereas 
the reference fiuid contains A^o particles and, in general, 

I3F,4p] = f3F,M - cl^\po) I dV[p(r)-po] 



dVdV'4')(|r-r'|;po) 



X [p(i-) -po][p(r') -Po] 



(24) 



(1'] f2) 

Using Cq (r;po) — Cq ' {r) = c(r) — —j3v{r) and the 
sum rule (fT4|) for n — 1 together with the vanishing of 
the excess chemical potential at zero density, we readily 
obtain 

c'i\po)^-f3mpo- (25) 

Formally substituting in Eq. ([M]) . po — > and p(r) po 
and making use of the fact that the excess free energy 
of a system vanishes with the density, we also obtain the 
dependence of the fluid excess free energy on the density: 



Nn 



(26) 



with the particle number Nq of the reference fluid and 
the Fourier transform v{k) of the interaction potential. 
Introducing Eqs. (|25p and ((26| into the Taylor expansion, 
Eq. ([M]) above, we obtain: 



+ 



^pmpo 

PviO)poiN-No 

Ppo 

^pmpo- 



dVdVw(|r-r'|)p(r)p(r') 
dVdV'y(|r-r'|)p(r) 



(27) 



Introducing x = r — r' the fourth term above becomes: 
-/3po / dVp(r) / d^xvil^l) = -PNpoiiO). (28) 
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Now the sum of the 1st, 2nd, 4th and 5th term in Eq. 
(PTj) yields: 



-PNpom 



^Po + Po{N- No) - Npo + ^po 



= Pv{0) [Nopo + Po{N- No) - Npo] - 0. (29) 

This is a remarkable cancelation because then only the 
3rd term in Eq. (|27p survives and we obtain: 



FeAp] 



d^rd^r'v{\r-r'\)p{r)p(r'), (30) 



which is our desired result.— i^^i^ 

The derivation above demonstrates that the excess free 
energy of any inhomogeneous phase for our ultrasoft flu- 
ids is given by Eq. ([30)l . irrespectively of the density of 
the reference fluid po- This is particularly important be- 
cause, usually, functional Taylor expansions are carried 
out around a reference fluid whose density lies close to 
the average density of the inhomogeneous system (crys- 
tal, in our case). However, in our systems this is impos- 
sible. The crystals occur predominantly in domains of 
the phase diagram in which the reference fluid is mean- 
ingless, because they are on the high density side of the 
A-line. It is therefore important to be able to justify the 
use of the functional and to avoid the inherent contradic- 
tion of expanding around an unstable fluid. In practice, 
of course, the higher-order dcf's do not exactly vanish, 
hence deviations from result are expected to occur, 
in particular at low temperatures and densities. Never- 
theless, the comparisons with simulations, e.g., in Refs. 
[l^l and [1^ fully justify our approximation. 

A mathematical proof of the mean-fleld character for 
fluids with inflnitely long-range and infinitesimally str ong 
repulsions has existed since the late 1970's, see Refs. [7y] 
and [77| . However, even far away from fulflUment of this 
limit, and for conditions that are quite realistic for soft 
matter systems, the mean-fleld behavior continues to be 
validi^i^^i^ The mean-fleld result of Eq. pO| has been 
put forward for the Gaussian model at high densitieSfSS, 
on the basis of physical argumentation: in the absence 
of diverging excluded-volume interactions, at sufflciently 
high densities any given particle sees an ocean of others - 
the classical mean-field picture. The mean-fleld character 
of the Gaussian model for moderate to high temperatures 
was demonstrated independently in Ref. [55]. Here, we 
have provided a more rigorous justification of its validity, 
based on the vanishing of high-order direct correlation 
functions in the fluid. It must also be noted that the 
mean-fleld approximation has recently been applied to a 
system with a broad shoulder and a much shorter hard- 
core interaction, providing good agreement with simula- 
tion results'^ - and allowing for the formulation of a gener- 
alized clustering criterion for the inhomogeneous phases. 



An astonishing similarity exists between the mean-fleld 
functional of Eq. (|30p and an exact result derived for 
inflnite-dimensional hard spheres. Indeed, for this case 
Frisch et a/| 78i79i80,8i .^gjj gag^jji and Rice^^ have 
shown that 



d^rd^r7(|r-r'|)p(r)p(r'), (31) 



where D ^ oo and /(|r — r'|) is the Mayer function of 
the inflnite dimensional hard spheres. Again, one has a 
bilinear excess functional whose integration kernel does 
not depend on the density; in this case, this is minus the 
(bounded) Mayer function whereas for mean-fleld fluids, 
it is the interaction potential itself, divided by the ther- 
mal energy ksT. In fact, the Mayer function and the di- 
rect correlation function coincide for inflnite-dimensional 
systems and higher-order contributions vanish there as 
weWr^ making the analogy with our three-dimensional, 
ultrasoft systems complete. Accordingly, inflnite dimen- 
sional hard spheres have an instability at some flnite k at 
the density p^, given by Pif{k) — 1. This so-called Kirk- 
wood instabilit'^ is of the same nature as our A-line but 
hard hyperspheres are athermal, so it occurs at a single 
point on the density axis and not a line on the density- 
temperature plane. Following Kirkwood's work;^ it was 
therefore argue d^°'^^ that hard hyperspheres might have 
a second-order freezing transition at the density ex- 
pressed as 

0.239(e/8)^/2gxp[^^(^/2)l/'^]i^l/^ 



P* 



(32) 

where the limit D oo must be taken and zo — 1.8558 is 
the value of the minimum of the Bessel function Jd/2{z) 
as -D — !■ oo. Note that p^, ^ as _D — s- oo. Later on, 
Frisch and Percus argued that most likely the Kirkwood 
instability is never encountered because it is preempted 
by a flrst-order freezing transitioni^ In what follows, we 
will show analytically that this is also the case for our 
systems, which might provide a flnite-dimensional real- 
ization of the above-mentioned mathematical limit. 



III. ANALYTICAL CALCULATION OF THE 
FREEZING PROPERTIES 

As mentioned above, an obvious candidate for a spa- 
tially modulated phase is a periodic crystal. The pur- 
pose of this section is to employ density functional the- 
ory in order to calculate the freezing properties. Under 
some weak, simplifying assumptions, the problem can be 
solved analytically. 

Adding the ideal contribution to the excess functional 
of Eq. ([50]) . the free energy of any spatially modulated 
phase is obtained as 



fceT J dV [In [p(r)A3] - l] 



dVdVw(|r-r'|)p(r)p(r'). (33) 
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As we are interested in crystalline phases, we parametrize 
the density profile as a sum of Gaussians centered around 
the lattice sites {R}, forming a Bravais lattice. In sharp 
contrast with systems interacting by hard, diverging po- 
tentials, however, the assumption of one particle per lat- 
tice site has to be dropped. Indeed, it will be seen that 
systems employ the strategy of optimizing their lat- 
tice constant by adjusting the number of particles per 
lattice site, Uc, at any given density p and temperature 
T. Accordingly, we normalize the profiles to ric and write 



R 



(34) 



R 



where the occupation variable Uc and the localization 
parameter a have to be determined variationally, and 
the lattice site density pi{r) is expressed as 



pi{r) = He (^^) 



3/2 



(35) 



Contrary to crystals of single occupancy, thus, the num- 
ber of particles N and the number of sites Ng of the 
Bravais lattice do not coincide. In particular, it holds 



N 



(36) 



and we are interested in multiple site occupancies, i.e., 
ric > 1 or even nc 3> 1; it will be shown that this cluster- 
ing scenario indeed minimizes the crystal's free energy. 

It is advantageous, at this point, to express the periodic 
density profile of Eq. ([34|) as a Fourier series, introducing 
the Fourier components of the same: 



P(r) 



K 



(37) 



where the set {K} contains all reciprocal lattice vectors 
(RLVs) of the Bravais lattice formed by the set {R}. 
Accordingly, the inverse of ([37]) reads asM- 



(38) 



PK = - dVe'K'-p(r) 

Vc Jc 



where the first integral extends over the elementary unit 
cell C of the crystal and Vc — V/Ng is the volume of 
C, containing a single lattice site. The second integral 
extends over all space, where use of the periodicity of p{r) 
and its expression as a sum over lattice site densities, Eq. 
([34ll . has been made. Using Eq. ((35|l we obtain 



PK 



(39) 



Note that the site occupancy Uc does not appear explic- 
itly in the functional form of the Fourier components of 
p{r), a feature that may seem paradoxical at first sight. 
However, for fixed density p and any crystal type, the lat- 
tice constant and thus also the reciprocal lattice vectors 
K are affected by the possibility of clustering, thus the 
dependence on Uc remains, albeit in an implicit fashion. 
With the density being expressed in reciprocal space, the 
excess free energy takes a simple form that reads as 



^ cx 



(40) 



K 



The ideal term, Fid [p] , can also be approximated ana- 
lytically, provided that the Gaussians centered at differ- 
ent lattice sites do not overlap. Let a denote the lattice 
constant of any particular crystal. Then, for aa^ 3> 1, 
the ideal free energy of the crystal takes the form 



N 



= In? 



3, faa' 
— m 

2 \ TT 



(41) 



where the trivial last term will be dropped in what fol- 
lows, since is also appears in the expression of the free 
energy of the fluid and does not affect any phase bound- 
aries. Putting together Eqs. and ([1T|) . we obtain 
a variational free energy per particle, /, for the crystal, 
that reads as 



Ne 



f{nc,a*-T\p*) 







5" 


Inric - 










" 2 



-yV(2"*) 



(42) 



where 



aa^ and Y = Kcr. In the list of arguments 



of / the first two are variational parameters whereas the 
last two denote simply its dependence on temperature 
and density. The free energy per particle, fsoi(T* , p*) = 
F/{Ne) of the crystal is obtained by minimization of /, 
i.e., 

/soi(T*,p*) = min T*,p*). (43) 

In carrying out the minimization, it proves useful to 
measure the localization length of the Gaussian profile, 
£ = in units of the lattice constant a instead of 

units of a. To perform this change, we first express the 
average density p of the crystal in terms of ric and a as 



P 



ZUc 
n3 



(44) 



where z is a lattice-dependent numerical coefficient of 
order unity. Introduce now the quantity aa^ — 7~^. 
Using Eq. (|^H) above, we obtain 



P 

ZTlc 



2/3 



(45) 
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This change of variables is just a mathematical trans- 
formation that simplifies the mathematics to follow; all 
results to be derived maintain their validity also in the 
original representation. For a further discussion of this 
point, see also Appendix B. 

Next we make the simplifying approximation to ig- 
nore in the sum over reciprocal lattice vectors on the 
right-hand-side of Eq. (|42p above all the RLVs beyond 
the first shell, whose length is Yi = Kia. This is justi- 
fied already because of the exponentially damping factors 
exp[— F^/(2a*)] in the sum. In addition, the coefficients 
|(5i(F)| themselves decay to zero as Y ^ oo, with an 
asymptotic behavior that depends on the form of 4>{x) in 
real space. The length of the first shell of RLVs of any 
Bravais lattice of lattice constant a scales as Ki = (/a, 
with some positive, lattice-dependent numerical constant 
C of order unity. Together with Eqs. this implies that 
the length of the first RLV depends on the aggregation 
number rir as 



ZUc 



1/3 



(46) 



and using Eq. (|45)) . we see that the ratio Yi /{2a*) takes 
a form that depends solely on the parameter 7, namely 



2a* 



2 



(47) 



Introducing Eqs. (gH]) and (gT]) into (g^]), we obtain 
another functional form for the variational free energy, 
/(ric, 7; T*, p*), expressed in the new variables. It can 
be seen that upon making the transformation (g5)l . the 
term 3/2 ln[(a*/7r)] delivers a contribution minus Inuc 
that exactly cancels the same term with a positive sign 
on the right-hand side of Eq. (gj). Accordingly, the only 
remaining quantity of the variational free energy that 
still depends on Uc is the length of the first nonvanishing 
RLV, Yi, whose ric dependence is expressed by Eq. (HS)) 
above. Putting everything together, we obtain 



fin,,r,T*,p*) 



In /?* — 1 — — [ln(77r) — 1] — In z 



^-lf~4>(Y,{n^))e-^^''\ 



(48) 



where is the coordination number of the reciprocal lat- 
tice. Minimizing / with respect to ric is trivial and using 
Eq. (gT]) we obtain 



5/ 

duc 



o^0'(yi)yi* = 0, 



(49) 



where the prime denotes the derivative with respect to 
the argument. Evidently, Yi coincides with the 
dimensionless wavenumber for which the dimensionless 
Fourier transform of the interaction potential attains its 



negative minimum. The other mathematical solution of 
(gni, Yi = 0, can be rejected because it yields nonpositive 
second derivatives or, on physical grounds, because it 
corresponds to a crystal with ric — > oo, whose occurrence 
would violate the thermodynamic stability of the system. 
Regarding second derivatives, it can be easily shown that 



and 



5V 



>0, 



Yi=v, 



djduc 



0, 



(50) 



(51) 



irrespective of 7. 

Having shown the coincidence of Yi with y^, we set 
4>{Yi) — (j){y^,) < in Eq. (|48|) above. Further, we notice 
that the term T* [In p* — 1] on the right-hand side of Eq. 

gives the ideal free energy of a uniform fiuid of den- 
sity p* and the term p* (}){{)) /2 the excess part of the same, 
see Eq. Subtracting, thus, the total fiuid free en- 

ergy pe_r particle, fn{T*,p*), we introduce the difference 
A/ = / — /fl, which reads as 



2tT* 

M{nc{y*),T,T*,p*)^ - — 



ln(77r) + 1 + 



21nz' 



^iP* It \ -7C'/2 

— -(/)(y*)e / . 



(52) 



The requirement of no overlap between Gaussians cen- 
tered on different lattice sites restricts 7 to be small; a 
very generous upper limit is 7 < 0.05. For such small val- 
ues of 7, the first term on the right-hand side of Eq. ((5^ 
above is positive. This positivity expresses the entropic 
cost of localization that a crystal pays, compared to the 
fiuid in which the delocalized particles possess transla- 
tional entropy. This cost must be compensated by a gain 
in the excess term, which is only possible if 0(j/*) < 0. An 
additional degree of freedom is offered by the candidate 
crystal structures. The excess free energy is minimized by 
the direct Bravais lattice whose reciprocal lattice has the 
maximum possible coordination number ^1. The most 
highly coordinated periodic arrangement of sites is fee, 
for which ^1 = 12. Therefore, in the framework of this 
approximation, the stable lattice is bcc. It must be em- 
phasized, though, that these results hold as long as only 
the first shell of RLVs is kept in the excess free energy. 
Inclusion of higher-order shells can, under suitable ther- 
modynamic conditions, stabilize fee in favor of bcc. We 
will return to this point later. 

Choosing now a as the edge-length of the conventional 
lattice cell of the bcc-lattice, we have z = 2 and C — 
2\/27r. Evidently, the lattice constant of the crystal is 
density- independent, a/a = (2-\/27r)/y,, contrary to the 
case of usual crystals, for which a oc The density- 

independence of a is achieved by the creation of clusters 
that consist of Uc particles, each of them occupying a 
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lattice site. The proportionality relation connecting ric 
and p* follows from Eq. (|^^ and reads as 



8\/27r3 

yl 



-p 



(53) 



It remains to minimize / (equivalently, A/) with re- 
spect to 7 to determine the free energy of the crystal. 
We are interested, in particular, in estimating the 'freez- 
ing line', determined by the equality of free energies of 
the fluid and the solid^i Accordingly, we search for the 
simultaneous solution of the equations 



dl 
A/ 



= 0, 
= 0, 



resulting into 



3r 



-T, — I — 0(y*)e 



and 



(t){.y*)e 



-7CV2 = 



BT* 



ln(77r) + 1 + 



21n; 



(54) 
(55) 

(56) 
(57) 



Substituting ([57)l into ((56|) and using z = 2 and C, — 
we obtain an implicit equation for 7 that reads as 



7-1 = -471^ 



In (77r) + 1 + 



2 In 2 



(58) 



and has two solutions, 71 = 0.018 and 72 = 0.038. Due to 
([50| and (I5ip , the sign of the determinant of the Hessian 
matrix at the extremum is set by the sign of f /d^^; 
Using Eqs. jM]), dST]), and ([58]) we obtain 



3T* _i 
972 27 



47r^ 



(59) 



which is positive for 7 = 71 but negative for 7 = 72. 
Only the first solution corresponds to a minimum and 
thus to freezing, whereas the second is a saddle point. 
Within the limits of the first-RLV-shell approximation, 
the crystals formed by Q^-potentials feature thus a uni- 
versal localization parameter at freezing: irrespective of 
the location on the freezing line and even of the interac- 
tion potential itself, the localization length ^ at freezing is 
a fixed fraction of the lattice constant and the parameter 
7 = {aa?)~^ attains along the entire crystallization line 
the value 



7f ^ 0.018. 



(60) 



We can understand the physics behind the constancy 
of the ratio l/a by examining anew the variational form 
of the free energy, Eq. ([42]) . Suppose we have a fixed den- 
sity p* and we vary ric, seeking to achieve a minimum of 
/. An increase in Uc implies an increase in the lattice con- 
stant a by virtue of Eq. The density profile takes 



advantage of the additional space created between neigh- 
boring sites and becomes more delocalized. This increase 
of the spreading of the profile brings with it an entropic 
gain which exactly compensates the corresponding loss 
from the accumulation of particles on a single site, ex- 
pressed by the term Inuc in Eq. (I42p . Expressing I in 
units of a, i.e., working with the variable 7 instead with 
the original one, a* , brings the additional advantage that 
7f becomes independent of the pair potential. The cor- 
responding value of a* at freezing, 0;^ , can be obtained 
from Eqs. pi)) and ([SH]). and reads for the bcc-lattice as 



87r27f 



(61) 



Here, a dependence on the pair interaction appears 
through the value of y*. 

Complementary to the localization parameter, we can 
consider the Lindemann ratio L at freezing)^^ taking into 
account that for the bcc lattice the nearest neighbor 
distance is d = a\/3/2. Employing the Gaussian den- 
sity parametrization, we find (r^) — 3/(2q:) and thus 
L = yj {r^)/d = y/2j. Using ((60| . the Lindemann ratio 
at freezing, Lf, is determined as 



Lf = 0.189. 



(62) 



This value is considerably larger than the typical value 
of 0.10 usually quoted for systems with harshly repulsive 
particles, such as, e.g., the bcc alcali metals and the fee 
metals Al, Cu, Ag, and Au but close to the value 
0.160 found by Stillinger and Weber— for the Gaussian 
core model. The particles in the cluster crystal are quite 
more delocalized than the ones for singly-occupied solids. 
The clustering strategy enhances the stability of the crys- 
tal with respect to oscillations about the equilibrium lat- 
tice positions. 

The locus of freezing points {T^,Pf) is easily obtained 
by Eqs. ([57)1 and ([55]) and takes the form of of a straight 
line: 



^ = 16^Sl^(y*)|e"'"'^'- 
Pi 



1.393 1(/)(2/,)|. (63) 



Contrary to the Lindemann ratio, which is independent 
of the pair potential, the freezing line does depend on the 
interaction potential between the particles. Yet, this de- 
pendence is a particular one, as it rests exclusively on the 
absolute value of the Fourier transform at the minimum, 
10(2/,)! and is simply proportional to it. Comparing with 
the location of the A-line from Eq. T^/pl = |0(y*)|, 
we find that crystallization preempts the occurrence of 
the instability: indeed, at fixed T*, pf < p\ or, equiva- 
lently, at fixed p* , > T{; see also Fig. Ill The transi- 
tion is first-order, as witnessed by the jumps of the values 
of a and pK at the transition, which are nonzero for the 
crystal but vanish in the fluid. This is analogous to the 
conjectured preemption of the Kirkwood instability for 
infinite-dimensional hard spheres by a first-order freez- 
ing transitioui^ 
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The freezing properties of Q^-potentials are, thus, 
quite unusual and at the same time quite simple: the 
lattice constant is fixed due to a clustering mechanism 
that drives the aggregation number ric proportional to 
the density. The constant of proportionality depends 
solely on the wavenumber for which the Fourier trans- 
form of the pair interaction has a negative minimum, Eq. 
(155)) . The freezing line is a straight line whose slope de- 
pends only on the value of the Fourier transform of the 
potential at the minimum, Eq. (j63p . The Lindemann 
ratio at freezing is a universal number, independent of 
interaction potential and thermodynamic state. 

Whereas the Lindemann ratio is employed as a mea- 
sure of the propensity of a crystal to melt, the height of 
the peak of the structure factor of the fluid is looked upon 
as a measure of the tendency of the fluid to crystallize. 
The Hansen- Verlet criterior&^ states that crystalliza- 
tion takes place when this quantity exceeds the value 
2.85. For the systems at hand, the maximum of S{y) lies 
at y,, as is clear from Eq. (PT|) . Using Eq. ([55]) for the 
location of the freezing line, we obtain the value 5'f (y,) 
on the freezing line as 



4>{y*) 



1.393 



3.542. 



(64) 



This value is considerably larger than the Hansen- Verlet 
threshold. In the fluid phase, Q^-systems can there- 
fore sustain a higher degree of spatial correlation before 
they crystallize than particles with diverging interactions 
do. This property lies in the fact that some contribu- 
tion to the peak height comes from correlations from 
within the clusters that form in the fluid; the formation 
of clusters already in the uniform phase is witnessed by 
the maxima of g{r) at r = seen in Fig. [1] and also 
explicitly visualized in our previous simulations of the 
modeli^ These, however, do not contribute to interclus- 
ter ordering that leads to crystallization. At any rate, 
the Hansen- Verlet peak height is also a universal quan- 
tity for all systems, in the framework of the current 
approximation. Moreover, both for the Lindemann and 
for the Hansen- Verlet criteria, the systems are more 
robust than usual ones, since they allow for stable fluids 
with peak heights that exceed 2.85 by 25% and for sta- 
ble crystals with Lindemann ratios that exceed 0.10 by 
almost 90%. 



COMPARISON WITH NUMERICAL 
MINIMIZATION 



IV. 



The density functional of Eq. (|33p is very ac- 
curate for the bounded ultrasoft potentials at 
hand,25^i55^^^^59^ rpj^g modeling of the in- 
homogeneous density as a sum of Gaussians is an 
approximation but, again, an accurate one, as has been 
shown by comparing with simulation resultsj^ see also 
section |V] of this work. The analytical results derived in 
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FIG. 4: The phase diagram of the GEM-4 model, obtained 
by full minimization of the density functional (|33p . under the 
Gaussian parametrization of the density, Eq. (|34p . redrawn 
from Ref. [2^. On the same plot, we show by the dotted line 
the approximate analytical result for the freezing line, Eq. 
(|63p . as well as the A- line of the system, Eq. (|22p . 



the preceding section rest on one additional approxima- 
tion, namely on ignoring the RLVs beyond the first shell. 
Here, we want to compare with a full minimization of 
the functional ([55)1 under the modeling of the density 
via (|34p , so as to test the accuracy of the hitherto drawn 
conclusions on clustering and crystallization. 

We work with the concrete GEM-4 system, for which 
the minimization of the density functional has been car- 
ried out and the phase diagram has been calculated in 
Ref. [1^. In Fig. U] we show the phase diagram ob- 
tained by the full minimization, compared with the freez- 
ing line from the analytical approximation, Eq. ()63|) . for 
this system. It can be seen that the latter is a very good 
approximation to the full result, its quality improving 
slowly as the temperature grows; the analytical approx- 
imation consistently overestimates the region of stabil- 
ity of the crystal. Moreover, whereas the approximation 
only predicts a stable bcc crystal, the high-density phase 
of the system is fee. Although bcc indeed is, above the 
triple temperature, the stable crystal immediately post- 
freezing, it is succeeded at higher densities by a fee lat- 
tice, which our analytical theory fails to predict. 

All these discrepancies can be easily understood by 
looking at the effects of ignoring the higher RLV shells 
from the summation in the excess free energy, Eq. (I42p . 
Consider first exclusively the bcc lattice. In Fig. [5] we 
show the locations of the bcc-RLVs, as obtained from the 
full minimization, by the downwards pointing arrows. It 
can be seen that the first shell is indeed located very 
closely to y,, as the analytical solution predicts. How- 
ever, the next two RLV shells do have contributions and, 
due to their location on the hump of (pijj), the latter is 
positive. By ignoring them in performing the analytical 
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FIG. 5: Inset: the Fourier transform (j){y) of the GEM-4 po- 
tential. Main plot: a zoom at the region of in which 
the first few nonvanishing RLV shells of the cluster-crystals 
of Fig.|4]lie. The arrows denote the positions of the shells and 
the numbers in square brackets the numbers of distinct RLVs 
within each shell. These positions are the result of the full 
minimization of the density functional. Downwards pointing 
arrows pertain to the direct bcc lattice and upwards pointing 
arrows to the direct fee one. In agreement with clustering pre- 
dictions, the positions of the RLVs are density-independent. 

solution, we are artificially lowering the free energy of the 
crystal, increasing thereby its domain of stability. 

The occurrence of a fcc-lattice that beats the bcc at 
high densities is only slightly more complicated to under- 
stand. A first remark is that the parameter a* increases 
proportionally to p* /T*, see the following section. Hence, 
the Gaussian factors from RLVs beyond the first shell, 
exp[—Yj^/{2a*)], i > 2, gain weight in the sum as den- 
sity grows. The cutoff for the RLV-sum is now provided 
rather by the short-range nature of (j){y) than by the ex- 
ponential factors. Due to the increased importance of the 
contributions from the i > 2-terms in the excess free en- 
ergy sum, the relative location of higher RLVs becomes 
crucial and can tip the balance in favor of fee, although 
the bcc-lattice has a higher number of RLVs in its first 
shell than the fee. In Fig. [5] we see that this is precisely 
what happens: the second RLV shell of the fee is located 
fairly close to the first. In the full minimization, both of 
them arrange their positions so as to lie close enough to 
?/*. Now, a total of fourteen 1st- and 2nd-shell RLVs of 
the fee can beat the twelve Ist-shell RLVs of the bcc and 
bring about a structural phase transformation from the 
latter to the former. 

The relative importance of the first and second neigh- 
bors is quantified by the ratio 

A = i^SM exp[-(y| - Y,')/{2a*)], (65) 
|0(yi)| ^ ' ^ 

where ^2 is the number of RLVs in the second shell. If 
(j){y) does not decay sufficiently fast to zero as y grows. 



then the fee lattice might even win over the bcc ev- 
erywhere, since then A could be considerable even for 
values of a* close to freezing, which are not terribly 
high. In fact, the penetrable sphere model (GEM-m with 
TO ~> 00) does not possess, on these grounds, a stable bcc 
phase at alli^ The prediction of the analytical theory on 
bcc stability has to be taken with care and is conditional 
to A being sufficiently small. A quantitative criterion 
on the smallness of A is model specific and cannot be 
given in general. The determination of the stable phases 
of the GEM-TO family and their dependence on to can be 
achieved by employing genetic algorithms'^ and will be 
presented elsewhere 

Notwithstanding the quantitative discrepancies be- 
tween the simplified, analytically tractable version of 
DFT and the full one, which are small in the first place, 
the central conclusion of the former remains intact: the 
RLVs of the crystals are density- independent. Whereas 
the analytical approximation predicts that the length of 
the first RLV shell coincides with , the numerical min- 
imization brings about small deviations from this predic- 
tion. However, by reading off the relevant values from 
Fig. [5l we obtain Yi = 5.625 for the bcc and Yi = 5.441 
for the fcc-lattice of the GEM-4 model. Comparing with 
the ideal value — 5.573, we find that the deviation 
between them is only a few percent. Clustering takes 
place, so that the lattice constants of both lattices remain 
fixed, a characteristic that was also explicitly confirmed 
by computer simulations of the modeli^ 



V. CONNECTION TO HARMONIC THEORY 
OF CRYSTALS 

The use of a Gaussian parametrization for the one- 
particle density profiles, Eq. ([34]), is a standard modeling 
of the latter for periodic crystals. This functional form is 
closely related to the harmonic theory of crystalsi^ Each 
particle performs oscillations around its lattice site, expe- 
riencing thereby an effective, one-particle site potential, 
Kiito(s) that is quadratic in the displacement s, for small 
values s/a [9l[. Here, we will explicitly demonstrate that 
the Gaussian form with the localization parameter pre- 
dicted from density functional theory coincides with the 
results obtained by performing a harmonic expansion of 
the said site potential. 

The formation of clustered crystals is a generic prop- 
erty of all Q^'^-systems, since the A-instability is com- 
mon to all of them; the form of the clusters that oc- 
cupy the lattice sites, however, can be quite complex, 
depending on the details of the interaction. The Gaus- 
sian parametrization p4p implies that for each of the ric 
particles of the cluster, the lattice site R is an equilib- 
rium position. In other words, the particular clusters we 
consider here are internally structureless. Clusters with 
a well-defined internal order have been found when an 
additional hard core of small extent is introducedp24 A 
necessary requirement for the lack of internal order is 
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that the Laplacian of the interaction potential v{r) be fi- 
nite at r = 0, as will be shown shortly. On these grounds, 
we impose from the outset on the interaction potential 
the additional requirement: 



oo for r ^ 0, (66) 



where the primes denote the derivative with respect to 
r. Eq. (|66p implies that v'{r) must be at least linear in r 
as r ^ 0. Concomitantly, v"{r) must be at least 0(1) as 
r ^ 0. As a consequence, we have 



v"ir) 



v'{r) 



for r — > 0. 



(67) 



It can be easily checked that (1571) is satisfied by all mem- 
bers of the GEM-m class for m > 2. It is also satisfied 
by the m = 2 member, i.e., the Gaussian model, which 
does not display clustering because it belongs to the (5+- 
class. This is, however, no contradiction. As mentioned 
above, the condition (j67p is necessary for the formation 
of structureless clusters and not a sufficient one. For Q^- 
potentials for which ([57]) is not fulfilled (such as the Fermi 
distribution models of Ref. [25j|), this does not mean that 
clusters do not form; it rather points to the fact that they 
possess some degree of internal order. 

The clustered crystals can be considered as Bravais 
lattices with a Uc-point basis. Accordingly, their phonon 
spectrum will feature 3 acoustic modes and 3(nc — 1) 
optical modes, for which the oscillation frequency u^k) 
remains finite as fc ^ 0. We are interested in the case 
Uc ^ 1, i.e., deep in the region of stability of the crystal, 
where the clusters have a very high occupation number. 
Consequently, the phonon spectra and the particle dis- 
placements will be dominated by the optical branches. 
Further, we simplify the problem by choosing, in the 
spirit of the Einstein model of the crystaljS one specific 
optical phonon with A: = as a representative for the 
whole spectrum. This mode corresponds to the relative 
partial displacement of two sublattices: one with Uc — 1 
particles on each site and one with just the remaining 
one particle per site. The two sublattices coincide at the 
equilibrium position and maintain their shape through- 
out the oscillation mode, consistent with the fact of an 
infinite- wavelength mode, k — 0. Accordingly, the site 
potential felt by any one of the particles of the single- 
occupied sublattice, K!ite(s), can be expressed as 



Kiitc(s) = (ric - 1) 



i;(s)+^t-(|s-R|) 

R#0 



(68) 



where s is the relative displacement of the two sublattices. 
For brevity, we also define 

_ "K3ite(s) 



W,ite(s) 



1 



(69) 



The Taylor expansion of a scalar function / around a 
reference point r reads as 

/(r + s) = /(r) + s • V/(r) + i (s • V)' /(r) + • • • (70) 



Setting r ^ and / Wsitc, we obtain the quadratic 
expansion of the site potential; the constant Vsitc(O) is 
unimportant. For the linear term, we have 



VM/sito(r-0) =w'(r)f 



R#0 



(71) 



where f and R are unit vectors. The sum in ([TTj) van- 
ishes due to lattice inversion symmetry; the first term 
also, since v'{0) = 0. Thus, the term linear in s in the 
expansion of \4ito(s) vanishes, consistently with the fact 
that s = is an equilibrium position. 

We now introduce Cartesian coordinates and write r = 
{x,y,z), s = {sx,Sy,Sz), and R = {R^, Ry, Rz)- The 
quadratic term in (|70p takes the explicit form 
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dxdy 



' dydz 
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dzdx 



2 (^'a^2 +^5^2 +^"5^2 ) • (72) 



Let us consider first the mixed derivative acting on 
Wsito(r), evaluated at r = 0. Using Eqs. and 
we obtain 







'xy / 


dxdy 


r=0 


r2 I, 
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R#0 



RxRy 

~R2- 



v'{r) 



v"{R) 



r=0 

v'{R) 



R 



.(73) 



The first term on the right-hand side of (|73p vanishes 
by virtue of (|^. The second one also vanishes due to 
the cubic symmetry of the lattice. Clearly, the other two 
terms in ()72|) with mixed derivatives vanish as well. For 
the remaining terms, the cubic symmetry of the lattice 
implies that all three second partial derivatives of Wsite 
at r = are equal: 



a2Wsite(r) 



i9a;2 



r=0 



9'W,ite(r) ^ d^W,^,{v) 

r=0 



1 



dy"^ 



V2w,ite(r) 



(}z2 



r=0 

(74) 



Gathering the results, we obtain the expansion of 
V^itc(s) to quadratic order in s as 



Fsitc(s) = V^sitc(0)-t- 



(»c - 1) 

6 



Ev^.(i?) 



R 



(75) 



which is isotropic in s, as should for a crystal of cu- 
bic symmetry. The one-particle motion is therefore har- 
monic; as we consider ^ 1, we set Uc — 1 = in (|75p 
and we obtain the effective, one-particle Hamiltonian Tii 
in the form: 



Hi 



P 

2m 



KS 



(76) 
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with 



(77) 



R 



and the momentum p and mass m of the particle. The 
density profile pi (r) of this single-particle problem is eas- 
ily calculated as {S{r — s))ni, yielding 



(78) 



This is indeed a Gaussian of a single particle, with a 
localization parameter ah = /3k; the total density on a 
given site will be then just ncPi(r), in agreement with 
the functional form put forward in Eq. (p4)) . 

It is useful to consider in detail the form of the local- 
ization parameter a^ predicted by the harmonic theory. 
The parameter k is expressed as a sum of the values of 
tpir) = V'^v{r) over the periodic set {R}. For every 

function ^/'(r) that possesses a Fourier transform ^p{k), it 
holds23 



(79) 



R 



K 



where {K} is the set of RLVs of {R} and ps = Ng/V is 
the density of lattice sites of {R}. From ([55]) . Ps"-c = P- 
Taking into account that the Fourier transform of V^w(r) 
is —kF'v(k), we obtain for the localization parameter of 
the harmonic theory the result 



ah 



(80) 



K 



The localization parameter must be, evidently, positive. 
Eq. (j80p manifests the impossibility for cluster formation 
if the Fourier transform of the pair potential is nonneg- 
ative, i.e., for interactions. In the preceding section, 
we showed within the DFT formalism that if the poten- 
tial is Q^, this implies the formation of cluster crystals. 
Harmonic theory allows us to make the opposite state- 
ment as well: if the potential is not Q^, then there can be 
no clustered crystals. Therefore, an equivalence between 
the character of the interaction and the formation 
of clustered crystals can be established. Moreover, Eq. 
([50)1 offers an additional indication as to why the RLV 
where v{K) is most negative is selected as the shortest 
nonvanishing one by the clustered crystals: this is the 
best strategy in order to keep the localization parameter 
positive. 

Harmonic theory provides, therefore, an insight into 
the necessity of locating the first shell of the RLVs at fc, 
from a different point of view than density functional the- 
ory does. The choice Ki = fc* guarantees that the parti- 
cles inhabiting neighboring clusters provide the restoring 
forces that push any given particle back towards its equi- 
librium position. The density functional treatment of the 
preceding sections establishes that the lattice constant is 
chosen by Q^-systems in such a way that the sum of 



intracluster and intercluster interactions, together with 
the entropic penalty for the aggregation of particles 
is optimized.— The unlimited growth of ric is avoided by 
the requirement of mechanical stability of the crystal. In- 
deed, for too high ric-values, the lattice constant would 
concomitantly grow, so that the resulting restoring forces 
working against the thermal fluctuations, would become 
too weak to sustain the particles at their equilibrium po- 
sitions. 

Let us, finally, compare the result ([80]) for the localiza- 
tion parameter with the prediction from DFT. We con- 
sider the high-density crystal phase, for which a* is very 
large, so that the simplification that only the first shell 
of RLVs can be kept in (|^ must be dropped. Setting 
df /da* = there, we obtain 



2a* {2a* f ^ ^ ' 



0. 



(81) 



The function (j){y) is short-ranged in reciprocal space, 
thus the sum in ([8T|) can be effectively truncated at some 
finite upper cutoff Vc- Then, there exists a sufficiently 
large density p* beyond which the parameter a* is so 
large that Y^/{2a*) < 1 for all Y < included in the 
summation. Accordingly, we can approximate all expo- 
nential factors with uniti in (I81|) . obtaining an algebraic 
equation for a*. Reverting back to dimensional quanti- 
ties, its solution reads as 



(82) 



K 



and is identical with the result from the harmonic theory, 
Eq. ([50]) . Thus, density functional theory and harmonic 
theory become identical to each other at the limit of high 
localization. This finding completes and generalizes the 
result of Archer,— who established a close relationship 
between the mean-field DFT and the Einstein model for 
the form of the variational free energy functional of the 
system. 

Consistently with our assumptions, a indeed grows 
with density. In fact, since the set of RLVs in the sum of 
(|82p is fixed, it can be seen that a is simply proportional 
to p/T. This peculiar feature of the class of systems we 
consider is not limited to the localization parameter, and 
its significance is discussed in the following section. 



VI. CONNECTION WITH INVERSE-POWER 
POTENTIALS 

As can be easily confirmed by the form of the den- 
sity functional of Eq. p3|l , the mean-field nature of the 
class of ultrasoft systems considered here (both Q^- and 
Q^-potentials) implies that the structure and thermody- 
namics of the systems is fully determined by the ratio 
p*/T* between density and temperature and not sepa- 
rately by p* and T*. This is a particular type of scaling 
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between the two relevant thermodynamic variables, rem- 
iniscent of the situation for systems interacting by means 
of inverse-power-law potentials v{r) having the form: 



v{r) 



'a 
. r 



(83) 



For such systems, it can be shown that their statistical 
mechanics is governed by a a single coupling constant 
T£>{n) expressed as2^ 



P 



(T*) 



(84) 



where D is the space dimension. It would appear that 
inverse-powers n = D satisfy precisely the same scaling 
as mean-field systems do but there is a condition to be 
fulfilled: inverse-power systems are stable against explo- 
sion, provided that n > D; this can be most easily seen by 
considering the expression for the excess internal energy 
per particle, u[p,T), given hy^ 



<P,T)^ 



pea- 



[{DI2)-1]\J, 



l{r-p,T)dr. (85) 



As g{r) 1 for r ^ oo, we see that the integral in ([SS]) 
converges only if n > D; a logarithmic divergence re- 
sults for n — D. A 'uniform neutralizing background' 
has to be formally introduced for n < Z?, to obtain 
stable pseudo one-component systems, such as the one- 
component plasmai^i^i^^ Since we aim at staying with 
genuine one-component systems throughout, we must 
strictly maintain n > D. 

Instead of taking the limit n — > D, we consider there- 
fore a different procedure by setting 



D + 6 



(86) 



with some arbitrary, finite S > 0. Then the coupling 
constant r£)(n) becomes 



Td{D + S)=p* (T 



(87) 



Now take Z? — > oo in this prescribed fashion, obtaining: 



lim Td(D 

D^oo 



S) 



£1 



which has precisely the same form as the coupling con- 
stant of our systems. In taking the limit D ^ oo, the 
exponent n = D + 6 of the inverse-power potential di- 
verges as well. It can be easily seen that in this case, the 
interaction v{r) of Eq. ([55)1 becomes a hard-sphere po- 
tential of diameter a. In other words, the procedure pre- 
scribed above brings us once more to infinite-dimensional 
hard spheres, because also the inverse-power interaction 
becomes arbitrarily steep. This is a v ery different way of 
taking the limit than in Refs. 80lf8lll83 |: there, the in- 
teraction is hard in the first place and subsequently the 
limit D — !■ cx) is taken, whereas here, interaction and di- 
mension of space change together, in a well-prescribed 
fashion. 



The fact that the statistical mechanics of ultrasoft flu- 
ids in three dimensions is determined by the same di- 
mensionless parameter as that of a particular realiza- 
tion of hard spheres in infinite dimensions is intrigu- 
ing. In a sense, ultrasoft systems are effectively high- 
dimensional, since they allow for extremely high densi- 
ties, for which every particle interacts with an exceed- 
ingly high number of neighbors. They might, in this 
sense, provide for three-dimensional approximate realiza- 
tions of infinite-dimensional models. This is yet another 
relation to infinite-dimensional systems, in addition to 
the one discussed at the end of Sec. |TT1 Whether there 
exists a deeper mathematical connection between the two 
classes, remains a problem for the future. 



VII. SUMMARY AND CONCLUDING 
REMARKS 

We have provided a detailed analysis of the proper- 
ties of bounded, ultrasoft systems, with emphasis on 
the Q^-class of interaction potentials. After having 
demonstrated the suppression of the contributions from 
the high-order direct correlation functions of the fluid 
phases (of order 3 and higher), we established as a 
consequence the accurate mean-fleld density functional 
for arbitrary inhomogeneous phases. Though this func- 
tional has been introduced and successfully used in the 
recent past both in statics^^i^SiSS^^^^^ g^^^ 

dynamicspi a sound justiflcation of its basis on the prop- 
erties of the uniform phase was still lacking. 

The persistence of a single, finite length scale for the 
lattice constants of the ensuing solids of Q^-systems has 
been understood by a detailed analysis of the structure of 
the free energy functional. In the fluid, the same length 
scale appears since the position of maximum of the liquid 
structure factor is independent of density. The negative 
minimum of the interaction potential in Fourier space 
sets this unique scale and forces in the crystal the forma- 
tion of clusters, whose population scales proportionally 
with density. The analytical solution of an approxima- 
tion of the density functional is checked to be accurate 
when confronted with the full numerical minimization 
of the latter. Universal Lindemann ratios and Hansen- 
Verlet values at crystallization are predicted to hold for 
all these systems, which differ substantially from those 
for hard matter systems. The analytical derivation of 
these results provides useful insight into the robustness 
of these structural values for an enormous variety of in- 
teractions. 

Though the assumption of bounded interactions has 
been made throughout, recent results^S indicate that 
both cluster formation and the persistence of the length 
scale survive when a short-range diverging core is super- 
imposed on the ultrasoft potential, provided the range 
of the hard core does not exceed, roughly, 20% of the 
overall interaction rangci^ The morphology of the result- 
ing clusters is more complex, as full overlaps are explic- 
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itly forbidden; even the macroscopic phases are affected, 
with crystals, lameUae, inverted lamellae and 'inverted 
crystals' showing up at increasing densities. The gener- 
alization of our density functional theory to such situa- 
tions and the modeling of the nontrivial, internal cluster 
morphology is a challenge for the future. Here, a mixed 
density functional, employing a hard-sphere and a mean- 
field part of the direct correlation function seems to be 
a promising way to proceediii^ Finally, the study of the 
vitrification, dynamical arrest and hopping processes in 
concentrated Q^-systems is another problem of current 
interest. The recent 'computer synthesis' of model, am- 
phiphilic dendrimers that do display precisely the form 
of Q^-interactions discussed in this work^ offers con- 
crete suggestions for the experimental realization of the 
hitherto theoretically predicted phenomena. 



APPENDIX B: PROOF OF THE EQUIVALENCE 
BETWEEN THE VARIATIONAL FREE 
ENERGIES / AND /. 

The introduction of the new variable 7 instead of a*, 
Eq. (|45|) . and the subsequent new form / of the varia- 
tional free energy, Eq. ([^5]) . are just a matter of conve- 
nience, which makes the minimization procedure more 
transparent. A free gift of the variable transformation 
is also the ensuing diagonal form of the Hessian matrix 
at the extremum. Fully equivalent results are obtained, 
of course, by working with the original variational free 
energy, /, Eq. (|42p . Here we explicitly demonstrate this 
equivalence. 

Keeping, consistently, only the first shell of RLVs with 
length Yi, / takes the form: 
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APPENDIX A: PROOF THAT THE GEM-m 
MODELS WITH m > 2 ARE Q±-POTENTIALS 

Consider the inverse Fourier transform of the spher- 
ically symmetric, bounded pair potential u(r), reading 
as 



v{r) 



1 

2^ 



k'v{k)- 



kr 



(Al) 



From (jAip , it is straightforward to show that the second 
derivative of v{r) at r = takes the form 



v"{r = 0) 



1 

6^ 



k%{k)dk. 



(A2) 



Evidently, if v"{r = 0) > 0, then v{k) must have negative 
parts and hence v{r) is Q^. For the GEM-m family, it 
is easy to show that v"{r = 0) = for m > 2, thus 
these members are indeed Q*, as stated in the main text. 
Double-Gaussian potentials of the form 
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where Yi and tIc are related via Eq. (|47p . Minimizations 
of / with respect to ric and a* yield, respectively: 



and 



-5^1/(2"*) 



0, (B2) 



T* + ^^(yi)r^2g-^iV(2«*) ^ 0. (B3) 

Subtracting the last two equations from one another we 
obtain 



^0' 
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(Yi)yie- 



0, 



(B4) 



implying Y\ = y*, as in the main text (once more, Yi = 
is a formal solution that must be rejected on the same 
grounds mentioned in the text.) From this property, Eq. 
([53|) immediately follows. Introducing A/ = / — /liq, we 
can determine a\ on the freezing line by requiring the 
simultaneous satisfaction of the minimization conditions 
above and of A/ = 0. The latter equation yields 





3 












- 1 






+2 









-^0(y.)e-^*/(^"^*), 



(B5) 



with |ei| > I £2 1, fi > (72, which feature a local minimum 
at r = are also Q^, for the same reason. Notice, how- 
ever, that v"{r = 0) > is a sufficient, not a necessary 



condition for membership in the 
exist Q^-potentials for which v" [r = 



-class. Thus, there 
0) < 0. 



which, together with Eq. (|B3p . yields, after some algebra 



In 



2/3' 



2a: 



(B6) 
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Using Eqs. (|45)) and (|46)) . we obtain the equation for the 
7f-parameter at freezing as 



hi (^z^/^'^itt'^ + 1 



(B7) 



which, upon setting z = 2 and C = 2\/2tt, yields Eq. 
of the main text. Alternatively, we can introduce the 
variable t = /y1 and rewrite Eq. (|B6p as 



In 



"^'^ Tit 



1 = 2t, 



(B8) 



which delivers t = 0.704 as a solution or, equivalently, 
t = (87r^7f)~^, in agreement with Eqs. ([60|) and ([6T|) of 
the main text. 
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